function th=theta(lat1,long1,lat2,long2)
    th=atan2(toRadians(lat2 - lat1), toRadians(long2 - long1) * cos(lat1))
end